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Abstract 

In a numerical Monte Carlo simulation of SU(2) Yang-Mills theory with 
light dynamical gluinos the low energy features of the dynamics as confine- 
ment and bound state mass spectrum are investigated. The motivation is 
supersymmetry at vanishing gluino mass. The performance of the applied 
two-step multi-bosonic dynamical fermion algorithm is discussed. 
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1 Introduction 



Supersymmctry seems to be a necessary ingredient of a quantum theory of gravity. It 
is generally assumed that the scale where supersymmctry becomes manifest is near to 
the presently explored clectrowcak scale and that the supersymmctry breaking is spon- 
taneous. An attractive possibility for spontaneous supersymmctry breaking is to ex- 
ploit non-perturbative mechanisms in supersymmetric gauge theories. Therefore the non- 
perturbative study of supersymmetric gauge theories is highly interesting [1]. In recent 
years there has been great progress in this field, in particular following the seminal papers 
of Seiberg and Witten [2] . 

The simplest supersymmetric gauge theories are the N — 1 supersymmetric Yang- 
Mills (SYM) theories. Besides the gauge fields they contain massless Majorana fermions 
in the adjoint representation, which are called gauginos in general. In the context of 
strong interactions one can call the gauge fields gluons and the gauginos gluinos. In the 
simple case of a gauge group SU(A'c) the adjoint representation is {N^ — 1) -dimensional, 
hence there are {N^ — 1) gluons and the same number of gluinos. 

The basic assumption about the non-perturbative dynamics of SYM theories is that 
there is confinement and spontaneous chiral symmetry breaking, similar to QCD. The 
confinement is realized by colourless bound states. Their mass spectrum is supposed to 
show a non-vanishing lower bound - the mass gap. Since external colour sources in the 
fundamental representation cannot be screened, the asymptotics of the static potential 
is characterized by a non-vanishing string tension. The expected pattern of spontaneous 
chiral symmetry breaking in SYM theories is quite peculiar: considering for definiteness 
the gauge group SU(A^c)) the expected symmetry breaking is Z2Nc ^2- For this we have 
recently found a first numerical evidence in a Monte Carlo simulation [3]. These general 
features of the low energy dynamics can be summarized in terms of low energy effective 
actions [4, 5]. 

The supersymmetric point in the parameter space corresponds to vanishing gluino 
mass {rrig = 0). For non-zero gluino mass the supersymmctry is softly broken and the 
physical quantities like masses, string tension etc. are supposed to be analytic functions of 
rrig. The linear terms of a Taylor expansion in rrig are often determined by the symmetries 
of the low energy effective actions [6]. 

The lattice regularization offers a unique possibility to confront the expected low 
energy dynamical features of supersymmetric gauge theories with numerical simulation 
results (for a recent review see [7].) On the lattice it is natural to extend the investigations 
to a general value of the gluino mass. In fact, to study exactly zero gluino mass is usually 
more difficult than the massive case and often an extrapolation to the massless point 
is necessary. The main difficulty in the numerical simulations is the inclusion of light 
dynamical gluinos. Although one can gain some insight also by studies of the quenched 
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theory [8, 9], the supersymmetry requires dynamical hght gluinos. 

In the present paper we report on a first large scale numerical investigation of SU(2) 
SYM theory with light gluinos. Although some preliminary results have already been 
published previously on different occasions [10, 11, 12, 13] and the question of the discrete 
chiral symmetry breaking has been dealt with in a recent letter [3] , this is the first detailed 
presentation of the obtained results. The numerical Monte Carlo simulations presented 
here have been performed on the CRAY-T3E computers at John von Neumann Institute 
for Computing (NIC), Jiilich. 

1.1 Lattice formulation 

For the lattice formulation we take the Wilson action, both for the gluon and gluino, as 
suggested some time ago by Curci and Veneziano [14]. The effective gauge field action is 

Scv^PT.(^-l^^Upi)-hogdetQ[U] . (1) 
pi 

For the gauge group SU(A^c), the bare gauge couphng is given by /3 = 2Nc/g^. The 
fermion matrix for the gluino Q is defined by 

4 

Qyv,xu = Qyv,xu[U] = ^yx^vu ~ ^ ['^2/,a;+/i ( 1 + 7//)^u,a;/i + ^y+il,xi}- ~ lfn)^vu,yiJ. ' (^) 

H=l 

K is the hopping parameter and the matrix for the gauge-field link in the adjoint repre- 
sentation is defined as 

Vrs,x, ^ Vrs,xAU] ^ 2T,{Ul%U,,%) = K.^., = K-^m " (3) 

The generators Tr = ^Xr satisfy the usual normahzation TV (A^As) = ^5rs- In case of 
SU(2) we have = with the isospin Pauli-matrices r^. 

The fermion matrix for the gluino Q in (2) is not hermitean but it satisfies 

= 75^75 ■ (4) 

This relation allows for the definition of the hermitean fermion matrix 

Q = ^,Q^qK (5) 

The factor | in front of logdet Q in (1) tells that we effectively have a flavour number 
Nf = I of adjoint fermions. This describes Majorana fermions in the Euclidean path inte- 
gral. For Majorana fermions the Grassmannian variables and are not independent 
but satisfy, with the charge-conjugation Dirac matrix C, 

= CC , = "^IC . (6) 
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(Note that here we use the Dirac-Majorana field instead of the Weyl-Majorana one 
A^.) The Grassmannian path integral for Majorana fermions is defined as 

J [d^]e-^*«* = y'[d^]e-^*^^«* = Pf (CQ) = Pf(M) . (7) 

Here the Pfaffian of the antisymmetric matrix M = CQ is introduced. The Pfaffian 
can be defined for a general complex antisymmetric matrix Map — —Mpa with an even 
number of dimensions (1 < a, /3 < 2A'") by a Grassmann integral as 

Pf (M) ^ I [dct>]e--^f"^-^t'^ = ^ea,p,...a.p,Ma,p, . . . M^^p, . (8) 

Here, of course, = d(l)2N ■ ■ ■ d4>i, and e is the totally antisymmetric unit tensor. It can 
be easily shown that 

[Pf (M)]^ = det M . (9) 

Besides the partition function in (7), expectation values for Majorana fermions can 
also be similarly defined [15, 9]. It is easy to show that the hcrmitean fermion matrix for 
the gluino Q has doubly degenerate real eigenvalues, therefore det Q = det Q = det M is 
positive and Pf(M) is real. In the effective gauge field action (1) the absolute value of 
the Pfaffian is taken into account. The omitted sign can be included by reweighting the 
expectation values according to 

- (signPf(M))cv ■ 

Here {■ ■ ■)cv means expectation value with respect to Scv- This sign problem is very 
similar to the one in QCD with an odd number of quark fiavours. 

The numerical simulations are almost always done on lattices with toroidal boundary 
conditions. In the three spatial directions it is preferable to take periodic boundary 
conditions both for the gauge field and the gluino. This implies that in the Hilbert space 
of states the supersymmetry is not broken by the boundary conditions. In the time 
direction in most cases we decided to choose periodic boundary conditions for bosons 
and antiperiodic ones for fermions, which is obtained if one writes traces in terms of 
Grassmann integrals: the minus sign for fermions is the usual one associated with closed 
fermion loops. 



1.2 Overview 

The aim of this paper is to give a complete presentation of the methods we used and 
to report on our numerical results. We concentrate on the confinement features as mass 
spectrum, emergence of supersymmetric multiplets of bound states and string tension. 

The plan of the paper is as follows: in the next section the numerical simulation algo- 
rithm is described. In particular, the computation of the necessary Pfaffians is dealt with 
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in section 2.3. The choice of algorithmic parameters and the observed autocorrelations 
are collected in section 2.5. The numerical results for the confinement potential and string 
tension as a function of the bare gluino mass are summarized in section 3. Section 4 is 
devoted to the spectrum of bound states: glueballs, gluino-glueballs and gluinoballs. This 
also includes the questions about possible mixing (section 4.4). The last section contains 
a summary. Three appendices are included: appendix A about least squares optimized 
polynomials, appendix B on the methods used for the determination of the smallest eigen- 
values and eigenvectors of the gluino matrix and appendix C about the main features of 
the C-|— I- implementation of our computer code. 



2 Multi-bosonic algorithm with corrections 

The multi-bosonic algorithm for Monte Carlo simulations of fermions has been proposed 
by Liischer [16]. In the original version for Nf (Dirac-) fermion flavours one considers the 
approximation of the fermion determinant 

|detW)r' = {detWtQ)}"'/^.g_2_, (U) 
where the polynomial P„ satisfies 

lim Pn{x) = x-^f/^ (12) 



n— »oo 



in an interval [e, A] covering the spectrum of Q^Q. Note that here the absolute value of 
the determinant is taken which leaves out its sign (or phase). In case oi Nf — ^, which 
corresponds to a Majorana fermion, this sign problem will be considered in section 2.3. 

For the multi-bosonic representation of the determinant one uses the roots of the 
polynomial rj, {j — 1, . . . ,n) 

n 

PniQ^Q) - PniQ') - ro UiQ' - r,) . (13) 

i=i 

Assuming that the roots occur in complex conjugate pairs, one can introduce the equiv- 
alent forms 

n n 
j=l j=l 

where rj = {fj^j + iuj)^ and pj = pj + iuj. With the help of complex scalar (pseudofermion) 
fields one can write 

ndet[(g-p*)(g-p,)]-i« m exp [(g-p*)(Q-p,)],.<|.,J . (15) 

j=l [ j=l xy J 

Since for a finite polynomial of order n the approximation in (12) is not exact, in 
principle, one has to extrapolate the results to n — > oo. In practice this can also be done 



by investigating the n-dependence and showing that the systematic errors introduced by 
the finiteness of n are neghgible compared to the statistical errors. 

The difficulty for small fermion masses in large physical volumes is that the condition 
number A/e becomes very large (10^ — 10^) and very high orders n — O(IO^) are needed 
for a good approximation. This requires large storage and the autocorrelation becomes 
bad since it is proportional to n. One can achieve substantial improvements on both these 
problems by introducing a two-step polynomial approximation [15, 17]. In this two-step 
multi-bosonic scheme (12) is replaced by 

4im^P£)(x)Pi^)(x) = , X e [6, A] . (16) 

The multi-bosonic representation is only used for the first polynomial P^^ which provides a 
first crude approximation and hence the order ni can remain relatively low. The correction 
factor P^^ is realized in a stochastic noisy correction step with a global accept-reject 
condition during the updating process (see section 2.1). In order to obtain an exact 
algorithm one has to consider in this case the limit n2 — > oo. For very small fermion (i.e. 
gluino) masses it turned out more practicable to fix some large n2 and perform another 
small correction in the evaluation of expectation values by reweighting with a still finer 
polynomial (see section 2.2). 



2.1 Update correction: global accept-reject 

The idea to use a stochastic correction step in the updating [18], instead of taking very 
large polynomial orders n, was proposed in the case of Nf — 2 fiavours in [19]. iV/ = 2 is 
special because the function to be approximated is just x'^ and P^f (x) can be replaced 
by the calculation of the inverse of xPj^^^ (x) . For general Nf one can take the two-step 
approximation scheme introduced in [15]. 

The two-step multi-bosonic algorithm is described in detail in [15]. Here we shortly 
repeat its main steps for the readers convenience and discuss the experience we obtained 
with it. The theory of the necessary optimized polynomials is summarized in appendix A 
following [17]. 

In the two-step approximation scheme for Nf fiavours of fermions the absolute value 
of the determinant is represented as 

|det(g)r^ - nr^^ ^r-^ • (17) 

detPiiHQ2)detPi^^(g2) 

The multi-bosonic updating with ni scalar pseudofermion fields is performed by heatbath 
and overrelaxation sweeps for the scalar fields and Metropolis sweeps for the gauge field. 
After a Metropolis sweep for the gauge field a global accept-reject step is introduced in 
order to reach the distribution of gauge field variables \U] corresponding to the right hand 
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side of (17). The idea of the noisy correction is to generate a random vector rj according 
to the normahzed Gaussian distribution 

e-'?^-Pn?(Q[t^]')'? 



(18) 



and to accept the change [W] <— [U] with probabihty 

min{l,A{rj;[U']^[U])} , (19) 

where 

Mv; [u'] - [u]) = exp {-v^Pil\Q[uf)v - v^Pil\Q[ur)v} . (20) 

The Gaussian noise vector rj can be obtained from r}' distributed according to the 
simple Gaussian distribution 

(21) 



J[dr]']e-^'^^' 
by setting it equal to 

V = Pi?iQ[U?)-^v' ■ (22) 
In order to obtain the inverse square root on the right hand side of (22), we can 
proceed with polynomial approximations in two different ways. The first possibility was 
proposed in [15] with x = as 

P(^)(x)-^ ^ ^ x'^f/'S^APiWx)] ■ (23) 

Here 

Sn.AP)^P'^ (24) 
is an approximation of the function on the interval P G [A^^^/^, e^^-'/^]. The poly- 
nomial approximations i?„3 and Sn^ can be determined by the same general procedure as 
P}^^ and Pj-^J . It turns out that these approximations are "easier" in the sense that for a 
given order higher precisions can be achieved than, say, for P^^^. 

Another possibility to obtain a suitable approximation for (22) is to use the second 
decomposition in (14) and define 

Pil^'^ {Q)^V^omQ- p.) . Pi? (Q') - PiT^ iQ)'PiT^ (Q) ■ (25) 

Using this form, the noise vector rj necessary in the noisy correction step can be generated 
from the gaussian vector rj' according to 

V = Pi^'J'KQrW : (26) 

where Pil^'^\Q)~^ can be obtained as 



= U2)Z[ ^ Pns{QWJ'\QV . (27) 
Pn2 [Q^) 
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In the last step P^g denotes a polynomial approximation for the inverse of Pj^ on the 
interval [e, A] . Note that this last approximation can also be replaced by an iterative 
inversion of PifiQ"^)- However, tests showed that the inversion by a least-squares opti- 
mized polynomial approximation is much faster because, for a given precision, less matrix 
multiplications have to be performed. 

In most of our Monte Carlo computations presented in this paper we used the second 
form in (26)-(27). The first form could, however, be used as well. In fact, for very high 
orders n2 or on a 32-bit computer the first scheme would be better from the point of 
view of rounding errors. The reason is that in the second scheme for the evaluation of 
PnJ'^\Q) we have to use the product form in terms of the roots pj in (25). Even using 
the optimized ordering of roots defined in [15, 17], this is numerically less stable than the 
recursive evaluation according to (63), (69). If one uses the first scheme both P^^ in (20) 
and i?„3 in (22)-(23) can be evaluated recursively. Nevertheless, on our 64-bit machine 
both methods worked well and we have chosen to apply (27) where the determination of 
the least-squares optimized polynomials is somewhat simpler. 

The global accept-reject step for the gauge field has been performed in our simulations 
after full sweeps over the gauge field links. The order rii of the first polynomial P^^J^ j^as 
been chosen such that the average acceptance probability of the noisy correction was 
near 90%. In principle one can decrease n\ and/or increase the acceptance probability 
by updating only some subsets of the links before the accept-reject step. This might be 
useful on lattices larger than our largest lattice 12^ • 24, but in our case we could proceed 
with full gauge sweeps and this seemed to be advantageous from the point of view of 
autocorrelations. 

2.2 Measurement correction: reweighting 

The multi-bosonic algorithms become exact only in the limit of infinitely high polynomial 

orders: n — > cxd in (12) or, in the two-step approximation scheme, 71,2 — cxd in (16). 
Instead of investigating the dependence on the polynomial order by performing several 
simulations, it is practically better to fix some high order for the simulation and perform 
another correction in the "measurement" of expectation values by still finer polynomials. 
This is done by reweighting the configurations in the measurement of different quantities. 
In case of Nj = 2 flavours this kind of reweighting has been used in [20] within the 
polynomial hybrid Monte Carlo scheme. As remarked above, Nj = 2 is special because 
the reweighting can be performed by an iterative inversion. The general case can, however, 
also be treated by a further polynomial approximation. 

The measurement correction for general Nf has been introduced in [12]. It is based 
on a polynomial approximation P^f which satisfies 

^lirn^ P£) (X)P(^) (a;)P£) (x) = x'^/Z^ , x e [e'. A] . (28) 
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The interval [e', A] can be chosen, for instance, such that e' = 0, A = Xmax, where Xmax is 
an absolute upper bound of the eigenvalues of Q^Q — Q^. In this case the limit ^4 — > oo 
is exact on an arbitrary gauge configuration. For the evaluation of P^^ one can use 714- 
independent recursive relations (see appendix A) , which can be stopped by observing the 
convergence of the result. After reweighting the expectation value of a quantity A is given 

by 

^ (Acxp{^t[i_p(4)(gtg)]^})^^^ 

where 77 is a simple Gaussian noise like r/' in (21). Here (. . .)^,, denotes an expectation 
value on the gauge field sequence, which is obtained in the two-step process described in 
the previous subsection, and on a sequence of independent t^'s. The expectation value 
with respect to the ?7-sequence can be considered as a Monte Carlo updating process with 
the trivial action = f]'^i]. The length of the ?7-sequence on a fixed gauge configuration 
can be, in principle, arbitrarily chosen. In praxis it has to be optimized for obtaining 
the smallest possible errors. If the second polynomial gives a good approximation the 
correction factors do not practically change the expectation values. A typical example is 
shown in figure 1. In such cases the measurement correction is good for the confirmation 
of the results. 

The application of the measurement correction is most important for quantities which 
are sensitive for small eigenvalues of the fermion matrix Q'^Q. The polynomial approxi- 
mations are worst near x — where the function x'^f^^ diverges. In the exact effective 
gauge action, including the fermion determinant, the configuration with small eigenvalues 
A are suppressed by K^sl'^. The polynomials at finite order are not able to provide such a 
strong suppression, therefore in the updating sequence of the gauge fields there are more 
configurations with small eigenvalues than needed. The exceptional configurations with 
exceptionally small eigenvalues have to be suppressed by the reweighting. This can be 
achieved by choosing e' = and a high enough order 714. It is also possible to take some 
non-zero e' and determine the eigenvalues below it exactly. Each eigenvalue A < e' is taken 
into account by an additional reweighting factor A.^f/'^P^^{K)P^^{A). The stochastic cor- 
rection in (29) is then restricted to the subspace orthogonal to these eigenvectors. Instead 
of e' > one can also keep e' = and project out a fixed number of smallest eigenvalues. 
Since the control of the smallest eigenvalues of the fermion matrix is an essential part of 
our simulations, a short summary of the numerical methods to obtain them is included 
in appendix B. 

Let us note that, in principle, it would be enough to perform just a single kind of 
correction. But to omit the reweighting does not pay because it is much more comfortable 
to investigate the (small) effects of different 77,4 values on the expectation values than to 
perform several simulations with increasing values of 722. Without the updating correction 
the whole correction could be done by reweighting in the measurements. However, in 
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practice this would not work either. The reason is that a first polynomial with relatively 
low order does not sufficiently suppress the exceptional configurations. As a consequence, 
the reweighting factors would become too small and would reduce the effective statistics 
considerably. In addition, the very small eigenvalues are changing slowly in the update 
and this would imply longer autocorrelations. 

A moderate surplus of gauge configurations with small eigenvalues may, however, be 
advantageous because it allows for a better sampling of such configurations and enhances 
the tunneling among sectors with different topological charges. For small fermion masses 
on large physical volumes this is expected to be more important than the prize one has to 
pay for it by reweighting, provided that the reweighting has only a moderate effect. The 
effect of a better sampling of configurations which small eigenvalues can be best illustrated 
by the distribution of quantities which diverge for zero eigenvalues. An example on 6^ • 12 
lattice at /5 = 2.3, K = 0.195 is shown in figure 2. 

2.3 The PfafRan and its sign 

The Pfaffian resulting from the Grassmannian path integrals for Majorana fermions (7) 
is an object similar to a determinant but less often used [21]. As shown by (8), Pf(M) is 
a polynomial of the matrix elements of the 2A'^-dimensional antisymmetric matrix M — 
—M^. Basic relations are 

M = P'^JP, Pf(M) = det(P) , (30) 

where J is a block-diagonal matrix containing on the diagonal 2(8)2 blocks equal to e = ia2 
and otherwise zeros. From this eq. (9) immediately follows. 

The form of M in (30) can be achieved by a procedure analogous to the Gram-Schmidt 
orthogonalization and, by construction, P is a triangular matrix. In order to see this, let 
us introduce the notation 

2N 

(uMv) = UaMapVp = (vM^u) (31) 
a,/3=l 

and denote the orthonormal basis vectors by {e^, a = 1,2,..., 2A^}. We are looking for 
a new basis {aj, bk, j, k = 1,2, . . . , N} obtained by 

% = P^2j-i = XI (ea-Pesj-i) , bk = Pe2k = X] (ca-Peafc) (32) 

a a 

such that the matrix elements on it are given by 

{ajMau) = , {bjMbk) = , (6feMa^) = -(ajMbk) = 5jk . (33) 
The construction is started by defining 

a^^e^, ^ ■ (34) 
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(If M21 is zero one has to rearrange the ordering of the original basis to achieve M21 7^ 0.) 
In the next step e;, Z = 3, 4, . . . , 2A'" is replaced by 

e'i_2 = ei-ai (hMci) + h (aiMei) , (35) 

which satisfy 

(el.Ma,) = (el^Mh) = . (36) 

With this the required form in (33) is achieved for oi and 61 and the corresponding matrix 
elements of P in (32), which are necessary for ai and bi, are determined. To proceed one 
has to return to (34) with {ca, a = 1, 2, ... , 2A^'} replaced by {e'„, a = 1, 2, ... , 2N — 2} 
and obtain the next (a, 6)-pair, until the whole space is exhausted. This gives a numerical 
procedure for the computation of P and the determinant of P gives, according to (30), 
the Pfaffian Pf(M). Since P is (lower-) triangular, the calculation of detP is, of course, 
trivial. 

This procedure can be used for a numerical determination of the Pfaffian on small 
lattices [12]. On lattices larger than, say, 4^ • 8 the computation becomes cumbersome due 
to the large storage requirements. This is because one has to store a full Q (8) Q matrix, 
with Q being the number of lattice points multiplied by the number of spinor-colour 
indices (equal to 4(A^^ — 1) for the adjoint representation of SU(A^c))- The difficulty of 
computation is similar to a computation of the determinant of Q with LC/-decomposition. 

Fortunately, in order to obtain the sign of the Pfaffian occurring in the measurement 
reweighting formula (10), one can proceed without a full calculation of the value of the 
Pfaffian. The method is to monitor the sign changes of Pf (M) as a function of the hopping 
parameter K. Since at ii' = we have Pf(M) = 1, the number of sign changes between 
K = and the actual value of K, where the dynamical fermion simulation is performed, 
determines the sign of Pf(M). The sign changes of Pf (M) can be determined by the flow 
of the eigenvalues of Q through zero. As remarked already in the discussion before (10), 
the fermion matrix for the gluino Q has doubly degenerate real eigenvalues therefore 

n/2 

detM = detg = n A- , (37) 
1=1 

where Aj denotes the eigenvalues of Q. This implies 

Q/2 _ n/2 _ 

|Pf(M)| = n |A^| , =^ Pf(M) = l[ ~K. (38) 
i=i 1=1 

The first equahty trivially follows from (9). The second one is the consequence of the 
fact that Pf(M) is a polynomial in K which cannot have discontinuities in any of its 
derivatives. Therefore if, as a function of an eigenvalue Aj (or any odd number of 
eigenvalues) changes sign the sign of Pf (M) has to change, too. We tested the sign of the 
Pfaffian in our Monte Carlo simulations by this spectral flow method. 
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As a representative example, let us consider the Monte Carlo runs on 6^ • 12 lattice 
for K — 0.19, 0.196, 0.20. The number of gauge configurations with negative Pfaffian 
in some representative subsets of the measured gauge configurations is given in table 1. 
The flow of the lowest eigenvalues with the hopping parameter is shown in some 
examples in figures 3, 4. The conclusion is that the probability of negative Pfaffians 
at most parameter values is negligible. Only at the largest hopping parameter, which 
corresponds to a negative gluino mass beyond the chiral phase transition [3], there is a 
somewhat larger fraction with negative Pfaffians but their effect on the averages is still 
smaller than the statistical errors. Therefore taking the absolute value of the Pfaffian, as 
in eq. (1), gives in the physically interesting points a very good approximation. 

Table 1: The fraction of Pfaffians with negative sign at (5 — 2.3 on 6^ • 12 lattice for 
different hopping parameters K. 



K 


7^ configs. 


#of P/(M)<0 


fraction 


0.19 


3840 (60x64) 





< 0.0003 


0.196 


5248 (82x64) 


14 


0.0027 


0.2 


2304 (36x64) 


69 


0.03 



2.4 Preconditioning 

The difficulty of numerical simulations increases with the condition number A/e charac- 
terizing the eigenvalue spectrum of fermion matrices on typical gauge field configurations. 
As it is well known, one can decrease the condition number by preconditioning. Even-odd 
preconditioning in multi-bosonic algorithms have been introduced in [22]. This turned 
out to be very useful in our simulations. 

For cvcn-odd preconditioning the hermitean fermion matrix Q is decomposed in sub- 
spaces containing the odd, respectively, even points of the lattice as 

A ^ / 75 -K-f^Moe \ , . 

For the fermion determinant we have 

det g = det Q , with g = 75 - K^-f^M^^M^^ . (40) 

The matrix Q'^ has a smaller condition number than Q^. 

The condition number and its fluctuations on different gauge configurations are do- 
minated by the minimal eigenvalue. An example of a comparison of the fluctuations of 
the lowest eigenvalue of and is shown in figure 5. As one sees, in this case the 
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mean of the smallest eigenvalue becomes about a factor 4 larger due to preconditioning. 
At the same time the largest eigenvalue becomes smaller, therefore the average condition 
number becomes about a factor 5 smaller. 

2.5 Parameter choice and autocorrelations 

The two-step multi-bosonic algorithm has several algorithmic parameters which can be 
tuned to achieve optimal performance. In fact, our experience shows that this tuning can 
bring substantial gain in efficiency. 

Polynomial degrees: In order to fix the polynomial degrees ni...4, in practice one 
performs trial runs using increasing values. At the same time, by observing the range of 
eigenvalues, one also obtains the interval [e. A] . The final value of rii is fixed by ensuring 
a high acceptance rate, around 90%, in the update correction step. n2 has to be large 
enough to keep the measurement correction small in important physical quantities. The 
final precision of the updating is set by n^, therefore the choice of rii and n2 does not 
infiuence the expectation values. For showing a typical example, in the upper part of 
figure 6 the polynomial approximation P^^^ and the product P^i^Pj^f are plotted in the 
interval [e. A]. The product Pil^PifPjif is displayed in the lower part of the figure. To 
the left (label a) the interval covers the range of the fiuctuating smallest eigenvalue, 
whereas to the right (label b) the function is shown in the range of fiuctuations of the 
last small eigenvalue which was determined explicitly. (In this case the correction factors 
were calculated from the eight smallest eigenvalues exactly and in the orthogonal subspace 
stochastically.) 

To fix the degree of the third polynomial, ns, we consider the probability pi of the 
system to jump between two identical configurations. In the limit ^ oo this prob- 
ability tends obviously to 1. In practice is increased till we get pi ~ 0.99, which is 
acceptable from the algorithm precision point of view, as one is convinced by comparing 
the expectation values. The choice of has to be tested, in principle, by observing its 
effect on the expectation values. Usually it is possible to choose already n2 so large that 
the measurement corrections with a substantially higher n4 are negligible compared with 
the statistical errors. 

The parameters of the numerical simulations at (3 = 2.3 are summarized in table 2. 
The runs with an asterisk had periodic boundary conditions for the gluino in the time 
direction T, the rest antiperiodic. K is the hopping parameter and [e, A] is the interval 
of approximation for the first three polynomials of orders ni 2,3, respectively. The fourth 
polynomial of order is defined on [0, A]. In the last two columns the number of per- 
formed updating cycles, respectively, the number of parallelly updated lattices (Niat) are 
given. 

Optimal ordering of the roots: The roots of P^J^ have to be always calculated. 
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Table 2: Parameters of the numerical simulations at (3 — 2.3. The notations are 
explained in the text. 



lattice: L,T 


K 


e 


A 


rii 


n2 


n?, 




updates 


Nlat 


6,12* 


0.16 


0.008 


3.2 


8 


32 


32 


- 


374400 


64 


6,12* 


0.17 


0.008 


3.2 


8 


32 


32 


- 


332800 


64 


6,12* 


0.18 


0.008 


3.2 


8 


32 


32 


- 


540800 


64 


6,12* 


0.185 


0.002 


3.4 


12 


32 


48 


- 


384000 


64 


6,12 


0.185 


0.002 


3.4 


16 


100 


150 


200 


550400 


64 


6,12* 


0.19 


0.0002 


3.5 


16 


60 


96 


- 


712800 


64 


6,12 


0.19 


0.0005 


3.6 


20 


112 


150 


400 


1487360 


64 


6,12* 


0.1925 


0.00003 


3.7 


22 


66 


102 


400 


1280000 


64 


6,12 


0.1925 


0.0001 


3.7 


22 


132 


180 


400 


3655680 


64 


6,12* 


0.195 


0.00003 


3.7 


22 


66 


102 


400 


1224000 


64 


D,iz 


u.iyo 


U.UUUUi 


6. i 




zUU 


onn 
oUU 


/inn 
4UU 


A cir\oc\c\ 
46UoU0 


64 


6,12 


0.196 


0.00001 


3.7 


24 


200 


300 


400 


952320 


64 


6,12 


0.1975 


0.000001 


3.8 


30 


300 


400 


500 


506880 


64 


6,12 


0.2 


0.000001 


3.9 


30 


300 


400 


500 


599040 


64 


8,16 


0.19 


0.00065 


3.55 


20 


82 


112 




1038400 


32 


8,16 


0.1925 


0.0001 


3.6 


22 


142 


190 




870400 


32 


12,24 


0.1925 


0.0003 


3.7 


32 


150 


220 


400 


216000 


9 



As discussed in section 2.1, depending on the way of doing the global accept-reject in 
updating, sometimes the roots of Pj^^ are also needed. Concerning this point a non- 
trivial question is how to order the roots when the representation (14) is used. Choosing 
this order naively leads to overflow and underflow problems because the product in (14) 
involves in general very different orders of magnitude. A good solution [17] is minimizing 
the maximal ratio of the values x"Pp(x) for x G [e, A], where Pp{x) denotes the partial 
product under consideration. This is in practice achieved by considering a discrete number 
of points in the interval {xi, . . . ,xjv} where N = 0{n). This gives in general sufficient 
numerical stability even for orders of many hundreds (see also the tests performed in [23]). 

Autocorrelations: During our simulations autocorrelations of different quantities 
were determined. Here we report on the analysis for the 12^ ■ 24 lattice at = 2.3, K = 
0.1925. Results for the 6'^ ■ 12 and the 8"^ ■ 16 lattices can be found in [3, 13]. We considered 
the short range exponential autocorrelation Texp of three different quantities, namely 

• the a-T]' propagator 

• the gluino-glue propagator 
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• the plaquette 

In the last case the data has been sufficient to give also an estimate of the integrated 
autocorrelation. 

Wc started the analysis by calculating the autocorrelation function C{t) for all these 
quantities. In case of the a-rj' we calculated the autocorrelation of the propagator at time 
distance At — 1, considering every 150'th configuration of the updating series. This was 
done separately on each of the lattices run in parallel. By averaging over the correlation 
functions obtained in this way we observed that the mean correlation function C°'~'^' [t] 
was qXj t — 1 already compatible with zero (C"-'' (1) = 0.028(19)), which lead to the 
conclusion that t^'^' < 150 updates. 

Estimating the exponential autocorrelation r^xp of the gluino-glue propagator we pro- 
ceeded similarly as in the a-r]' case. On all nine lattices that were run in parallel we 
determined independently the autocorrelation of the propagator at time distance At — 1 
on every 150'th configuration of our total history. The exponential autocorrelation time 
was then estimated by fitting an exponential of the form ex.^{—t / Texp) to the first points 
of the curve for each lattice. A typical autocorrelation function with the exponential 
fit can be seen in figure 7. By finally taking the average of Texp over all lattices we ar- 
rived at the result displayed in table 3. It has to be understood that Texp determined in 
this way displays a mode between the true short range exponential and the integrated 
autocorrelation, since only every 150'th sweep has been considered. 

To estimate the integrated autocorrelation time Tint of the plaquette we proceeded in 
a different manner. On the basis of prior analysis [3, 13] we expect the order of magnitude 
of Tint to be 10^ ~ 10^. Since for each lattice we have a total of about 24000 configurations 
in equilibrium we expect our time history to have a length of at most ~ 100rj„(. This leads 
to the conclusion that standard methods to determine Tint [24, 25] are not reliable since 
they require statistics that are at least of the order of several hundred Tint- Therefore, 
to estimate the order of magnitude of Tint we proceed as follows. For each lattice run 
in parallel wc calculated the autocorrelation function C^''"''^{t) of the plaquette for the 
complete history of 24000 configurations. We fitted an exponential decay to C^^°''^{t) in 
a small interval (typically [0,300]) at the beginning where the fastest decay mode should 
be dominant. For longer distance the exponential typically decayed faster than C'^''"''^ (t) . 
This expected behaviour could usually be observed up to a point i where C^^'^'^it) started 
to be dominated by its noise. We then calculated 

t 

J2CP^'^{t) (41) 
t=l 

and took this value as an estimate of Tint- We expect this procedure only to lead to an 
order of magnitude estimate for the integrated autocorrelation. The typical behaviour 
for the autocorrelation function of the plaquette together with the exponential fit can be 
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observed in figure 8. In this example the cutoff i has been chosen at about i ~ 3500 
updates since at this point C^'"^(i) is clearly dominated by its noise. The final result for 
Texp and Tint have been obtained by averaging over all nine lattices run in parallel, and 
can be found in table 3. 

Table 3: Autocorrelation and integrated autocorrelation of the propagators and the 
plaquette on 12^ ■ 24 lattice at 13 = 2.3 and K = 0.1925. 







' ml 


a-Tj' 


< 150 




gluino-glue 


620(60) 


1100(200) 


plaquette 


378(37) 


675(43) 



3 Confinement potential 

The potential between static colour sources in gauge field theory is a physically very 
interesting quantity because it is characteristic for the dynamics of the gauge fields. If 
the sources are in the fundamental representation of the gauge group they can be called 

static quarks. 

For a model containing dynamical matter fields in the fundamental representation, 
as is the case for QCD with dynamical quarks, they will screen the static quarks. The 
potential then approaches a constant at large distances [26]. The string tension a, which 
is the asymptotic slope of the potential for large distances, vanishes accordingly. This 
type of screening is of a more kinematical nature. 

On the other hand, if only matter fields in the adjoint representation of the gauge 
group are present, as in the case of supersymmetric N=l Yang-Mills theory, there are 
different possibilities. Either the string tension does not vanish and static quarks are 
confined, or the static quarks are screened dynamically by the gauge fields. The latter sit- 
uation is found in two-dimensional supersymmetric Yang-Mills theory [27]. The screening 
mechanism is related to the chiral anomaly and appears to be specific to two dimensions. 

Four- dimensional SUSY Yang-Mills theory is believed to confine static quarks [28]. 
Furthermore, the behaviour of the string tension as a function of the gluino mass can give 
indications on the question, whether QCD and Super-QCD are smoothly hnked [29]. 

We have determined the static quark potential and the string tension for N=l SUSY 
Yang-Mills theory from our Monte Carlo results. The starting point are expectation 
values of rectangular Wilson loops {W{R,T)). In order to improve the overlap with 
the relevant ground state we have apphed APE-smearing [30] to the Wilson loops. The 
optimal smearing radius turns out to be near Rg — 3. 
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Prom the Wilson loops the potential can be found via 



V{R)=lim V{R,T) , (42) 

where 

V{R, T) = \og{W{R, T)) - \og{W{R, T + l)) . (43) 

The large T limit is approached exponentially [31]. We have obtained the potential V{R) 
through a fit of the form 

V{R, T) = V{R) + ci(i?)e-^^^^^^. (44) 

As an example we show \^(3,T) as a function of T on a 12^ ■ 24 lattice in figure 9. For 
T > 6 the errors grow significantly and we have chosen 1 < T < 6 as the best fit interval 
on this lattice. On the 8^ • 16 lattice fit intervals from T = 1 to 4 or 5 yield consistent 
results. 

In this way the static potential V{R) has been obtained for 1 < < 6 on the 8'^ ■ 16 
lattice and for 1 < i? < 9 on the 12'^ ■ 24 lattice. For larger values of R the errors become 
rather large and the results arc not reliable anymore. Anyhow, for R > L/2 increasing 
finite size effects arc to be expected. In figures 10 and 11 the potential is shown on the 
L = 8 lattice aX K = 0.19 and L = 12 lattice at K = 0.1925, respectively. 

The string tension a is finally obtained by fitting the potential according to 

(y 

V{R)^Vo-- + aR. (45) 
rC 

The value of a depends on the range of R taken for the fit. In general it tends to decrease 
if the largest values of R are included in the fit. However, this should not be interpreted 
as a signal for screening, since the potential is expected to bend down due to finite size 
effects. In table 4 the values for ^/a in lattice units are shown for different fit ranges. 



Table 4: Square root of the string tension a in lattice units and Coulomb strength a 
from fits to V{R) = ~ % + crR over different ranges of R. 



lattice 


K 


R fit range 


a^Ja 


a 


8^ • 16 


0.19 


1-4 


0.22(1) 


0.23(2) 


8^ • 16 


0.19 


1 - 5 


0.21(1) 


0.25(1) 


8^ • 16 


0.1925 


1-4 


0.21(1) 


0.23(2) 


8^ • 16 


0.1925 


1 - 5 


0.19(1) 


0.25(2) 


12^ . 24 


0.1925 


1 - 6 


0.17(1) 


0.25(2) 


12^ • 24 


0.1925 


1 - 7 


0.16(1) 


0.26(2) 


12 -^ • 2 1 


0.1925 


1 8 


0.13(2) 


0.31(1) 
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We consider the range l<i?<L/2as reliable and quote as final results for the string 
tension 



ax/a = 0.22(1) for = 0.1900, L = 8, 
= 0.21(1) for = 0.1925, L = 8, 
av^ = 0.17(1) for K = 0.1925, L = 12. 



(46) 



The string tension in lattice units is decreasing when the critical line is approached, as 
it should be. This is mainly caused by the renormalization of the gauge coupling due to 
virtual gluino loop effects which are manifested by decreasing lattice spacing a. Prom a 
comparison of the L — % and L = 12 results one sees that finite size effects still appear 
to be sizable. This has to be expected because we have for the spatial lattice extension 
L — 12a the result Ly^ ~ 2.1. In QCD with A^a ~ 0A5GeV this would correspond to 
L ~ 1 fm. Although we are deahng with a different theory where finite size effects as a 
function of L^/a are different, for a first orientation this estimate should be good enough. 

The coefficient a of the Coulomb term is close to the universal Liischer value of 
7r/12 = 0.26 [32]. 

For the ratio of the scalar glueball mass m(0+), to be discussed below, and the square 
root of the string tension we get 



The uncertainties are not very small, but the numbers are consistent with a constant 
independent of K in this range. They are of the same order of magnitude but somewhat 
smaller than in pure SU(2) gauge theory [33], where at /? = 2.5 we have m{0^)/y/a = 3.6- 
3.8, depending on the lattice size. 

4 Light bound state masses 

The non- vanishing string tension observed in the previous section is in accordance with the 
general expectation [1, 4] that the Yang-Mills theory with gluinos is confining. Therefore 
the asymptotic states are colour singlets, similarly to hadrons in QCD. The structure of 
the light hadron spectrum is closest to the (theoretical) case of QCD with a single fiavour 
of quarks where the chiral symmetry is broken by the anomaly. 

Since both gluons and gluinos transform according to the adjoint (here triplet) rep- 
resentation of the colour group, one can construct colour singlet interpolating fields from 
any number of gluons and gluinos if their total number is at least two. Experience in 
QCD suggests that the hghtest states can be well represented by interpolating fields build 



m(0+)/v^ 



3.4(7) 
3.0(4) 
3.1(7) 



for K = 0.1900, L = 8, 
for K = 0.1925, L = 8, 
for K = 0.1925, L = 12. 



(47) 
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out of a small number of constituents. Simple examples are the glueballs known from 
pure Yang-Mills theory and gluinoballs corresponding to pseudoscalar mesons. We shall 
call the simplest pseudoscalar gluinoball made out of two gluinos the a-rj' state. Here 
the label a reminds us to the fact that the constituents are in the adjoint representation 
and 77' stands for the corresponding Ty'-meson in QCD. Mixed gluino-glueball states can be 
composed of any number of gluons and any number of gluinos, in the simplest case just 
one of both. 

In general, one has to keep in mind that the classification of states by some interpo- 
lating fields has only a limited validity, because this is a strongly interacting theory where 
many interpolating fields can have important projections on the same state. Taking just 
the simplest colour singlets can, however, give a good qualitative description. 

In the supersymmetric limit at zero gluino mass rrig = the hadronic states should oc- 
cur in supermultiplets. This restricts the choice of simple interpolating field combinations 
and leads to low energy effective actions in terms of them [4, 5] . For non-zero gluino mass 
the supersymmetry is softly broken and the hadron masses are supposed to be analytic 
functions of rrig. The linear terms of a Taylor expansion in rrig are often determined by 
the symmetries of the low energy effective actions [6]. 

The effective action of Veneziano and Yankielowicz [4] describes a chiral supermultiplet 
consisting of the 0~ gluinoball a-rj', the 0"*" gluinoball a-/o, and a spin | gluino-glueball. 
There is, however, no a priori reason to assume that glueball states are heavier than 
the members of the supermultiplet above. Therefore Farrar, Gabadadze and Schwetz 
[5] proposed an effective action which includes an additional chiral supermultiplet. This 
multiplet consists of a 0+ glueball, a 0~ glueball and another gluino-glueball. The effective 
action allows mass mixing between the members of the two supermultiplets. The masses 
of the lightest bound states and the mixing among them can be investigated by Monte 
Carlo simulations. 

4.1 Glueballs 

The glueball states as well as the methods to compute their masses in numerical Monte 
Carlo simulations are well known from pure gauge theory. (For a recent summary of 
results and references see [33].) 

The lightest state is the = 0^ glueball which can be generated by the symmetric 
combination of space-like plaqucttes touching a lattice point. In order to optimize the 
signal and enhance the weight of the lightest state one is taking blocked [34] or smeared 
[30] links instead of the original ones. In order to obtain the masses, for a first orientation, 
one can use effective masses m(ti, ^2, T) assuming the dominance of a single state for time- 
slices ti, t2 on the periodic lattice with time extension T. One can search for time distance 
intervals where the effective masses are roughly constant and then try single mass fits in 
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these intervals. In cases with high enough statistics and corresponding small statistical 
errors two-mass fits in larger intervals can also be stable and give information on the mass 
of the next excited state. 

Since no previous results on the glueball mass spectrum with dynamical gluinos are 
available in the literature, we started our search for dynamical gluino effects on small 
lattices as 4^ • 8 at hopping parameter values K > 0.16. Wc observed some effects for 
K > 0.18 where we started runs on larger lattices, up to 12^ • 24. As already seen in 
the previous section, the lattice spacing a is decreasing with increasing K (i.e. decreasing 
gluino mass). This means that effectively we are closer to the continuum limit at larger 
K, resulting in smaller glueball (and other) masses in lattice units. This effect is strongest 
at zero gluino mass where a first order phase transition is expected due to the discrete 
chiral symmetry breaking. First numerical evidence for this phase transition has been 
reported by our collaboration at K = Kq = 0.1955(5) [3]. 

With our spectrum calculations wc stayed below this value and stopped at = 0.1925 
where the 12'^ -24 lattice is already not very large. The obtained masses for the 0"*" glueball 
in lattice units are 



am(0+) 


= 0.95(10) 


for K 


= 0.1800, L = 


6, 


am(0+) 


= 0.85(6) 


for K 


= 0.1850, L = 


6, 


am(0+) 


= 0.75(6) 


for K 


= 0.1900, L = 


8, 


am{0'^) 


= 0.63(5) 


for K 


= 0.1925, L = 


8, 


am(0"'") 


= 0.53(10) 


for K 


= 0.1925, L = 


12 



In addition to the — 0+ glueball we have studied the pseudoscalar 0~ glueball. In 
order to create a pseudoscalar glueball from the vacuum with an operator built from closed 
loops on the lattice, one needs loops which cannot be rotated into their mirror images. 
For gauge group SU(2) the traces of loop variables are real and do not distinguish the 
two orientations of loops. The smallest loops with the desired property are made of eight 
links. One possibility would be to take the simplest lattice version of Tr:{enj^pcrF'^^ Ff"^). 
However, it contains two orthogonal plaquettes and cannot be put into a single time-slice. 
Therefore we have chosen to take the loop C shown in figure 12 [35]. 

The time-slice operator for the pseudoscalar glueball is then given by 

S{t) = Y.[TiU{C)-TiU{PC)] , (49) 

R 

where the sum is over all rotations it! in the cubic lattice group and PC is the mirror 
image of C. As usual, APE-smearing has been applied to the links appearing in the loop. 

The pseudoscalar glueball mass has been calculated from the time-slice correlation 
functions as an effective mass from distances 1 and 2 with optimized smearing radius. On 
the 6^-12 lattice a good smearing radius is obtained for i?s = 4 or 5, and the numbers are 
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very stable. On the 8^-16 lattice a clear plateau in the number of smearing steps could 
not be seen. Nevertheless, for a smearing radius between 5 and 8 we obtain rather stable 
results. The masses in lattice units are 



am(0~) 


= 1.5(3) 


for K = 


0.1850, L = 


6, 


am(0"") 


= 1.45(10) 


for K = 


0.1900, L = 


6, 


am(0~) 


= 1.3(1) 


ior K ^ 


0.1925, L = 


6, 


am(0~) 


= 1.1(1) 


for X = 


0.1925, L = 


8. 



The pscudoscalar glueball appears to be roughly twice as heavy as the scalar one. This 
is similar to pure SU(2) gauge theory, where m(0~)/m(0"'") = 1.8(2) [33]. 



4.2 Gluino-glueballs 

One can construct colour singlet states from the gluinos and the field strength tensor in 
the adjoint representation. One of these states is a spin | Majorana fermion which occurs 
in the construction of the Veneziano- Yankielowicz effective action [4] . In order to find the 
lowest mass in this channel we consider the correlator consisting of plaquettes connected 
by a quark propagator line: 

^ggip^-iU) ~ "^"^s^yinov {XxQxr,ysXy) (51) 

where 

X^ = ^'IVcolour(T.C7.) (52) 

and the plaquette variable is defined as 

C7, = U^i{x, 12) + U^i{x, 13) + U^i{x, 23) . (53) 

For antiperiodic boundary conditions for the gluino in the time direction the correlator is 
antiperiodic. By inserting 74 into the correlation function (51) it becomes periodic also 
with antiperiodic boundary conditions. The resulting projection on the ground state have 
in both cases either been compatible with one another, or the propagator modified with 
74 has shown more mixing with larger masses. Therefore in extracting the masses we 
considered only the above propagator without 74. 

For the gluino-glueball, in order to obtain a satisfactory signal, APE-smearing [30] 
has been implemented for the hnks and Jacobi-smearing [36] for the gluino field. Tests 
have shown [13] that Teper-blocking for the hnks was in this case not as well suited. 
Table 5 shows the smearing parameters used for the gluino-glueball on different lattices 
at different hopping parameters. They have been optimized by measuring the masses on 
a small sample of data and tuning the parameters accordingly to obtain the lowest mass 
values. 
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Table 5: Smearing parameters for Jacobi and APE-smearing used for measuring the 
gluino-gluehall. 



lattice 


K 


■^Jacobi 


K Jacobi 


Nape 


^APE 


8=^-16 


0.19 


20 


0.22 


8 


0.35 


8=^ • 16 


0.1925 


23 


0.185 


12 


0.34 


12^ ■ 24 


0.1925 


19 


0.20 


9 


0.3 



The masses for the gluino-glueball were determined first by considering effective 
masses m(ti,t2,T) assuming the dominance of a single state for time-slices ti,t2 on the 
periodic lattice with time extension T. From this time distance intervals were determined 
where the effective masses were roughly constant and single mass fits in these intervals 
were performed. The results are shown in table 6. 

Table 6: Lowest masses for the gluino-glueball at different hopping parameters and 
lattices. The value of the gauge coupling has been {3 = 22 throughout. 



gluino-glueball 


X- 0.18 


X = 0.185 


X-0.19 


X = 0.1925 


6^.12 


1.93(5) 


1.39(8) 


1.05(20) 




8=^-16 






0.87(13) 


0.82(18) 


12=^ • 24 








0.93(8) 



4.3 Gluinoballs 

Besides the gluino-glueball in this work we consider also gluinoballs defined by a colourless 
combination of two gluino fields. The a-rj' has spin-parity 0^ and the a-/o spin-parity O"*". 
In the simulations for a-rj' and a-/o, respectively, the wave functions ^75 \E' and 'l'^ were 
used. These gluinoballs are contained in the Veneziano-Yankielowicz super- multiplet [4]. 
For the correlation function a straightforward calculation as in [15] with F G {1, 75} yields 

T-,-,{x, y) = (IV,, {rQ-i}iv,, {rg-i} - 2 iv,, {rg-^FQ-;}) . (54) 

Note the factor of two originating from the Majorana character of the gluinos. In ana- 
logy with a fiavour singlet meson in QCD the propagator consists of a connected and a 
disconnected part: the left, respectively, the right term of (54). 

The numerical evaluation of the time-slice of the connected part can be reduced to the 
calculation of the propagator from a few initial points. The disconnected part is calculated 
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using the volume source technique [37] . For the determination of the gluinoball propagator 
no smearing has been used. 

In case of the a-fo particle the disconnected and the connected parts are of the same 
order of magnitude. The former has a much worse signal to noise ratio than the latter. 
This leads to a larger error on the a-fo as compared to the a-rj' which is dominated by 
the connected part. 

Our results for the a-rj' and the a-fo masses for different lattices and hopping pa- 
rameters can be found in table 7. In case of the a-r}' the data has been good enough to 

Table 7: Lowest masses for the a-r]' and the a-fo at different hopping parameters 
and lattices. The gauge coupling is given by f] = 2.3 throughout. In the last column 
with a star the next higher mass is shown, whenever it could be determined. 



a-rj' 


X-0.18 


0.185 


X = 0.19 


K = 0.1925 


K = 0.1925* 




1.155(11) 


0.941(8) 


0.594(14) 






8=^-16 






0.725(20) 


0.551(17) 


1.282(26) 


12^ • 24 








0.48(5) 


1.09(5) 


a-fo 


K = 0.18 


K = 0.185 


K = 0.19 


K = 0.1925 




6^ ■ 12 


1.49(13) 


1.11(17) 








8^ ■ 16 






1.20(22) 


0.81(17) 




12^ • 24 








1.00(13) 





estimate also the next higher state. (These data can be found in the column denoted by 
a star.) The lowest masses have been obtained by using effective masses and fits as for 
the gluino-glueball. The fits were rather stable in case of the a-rj' on the 12^ • 24 lattice. 
This allowed to extract two masses from the data. Errors were estimated by the jackknife 
method. 

4.4 Glueball-gluinoball mixing 

In the low energy effective action of Farrar, Gabadadze and Schwetz [5] there is a possible 
non-zero mixing between the states in the two light supermultiplets. In particular there 
can be mixing of the a-fo gluinoball and the 0+ glueball. 

In order to study the mixing we have calculated the connected cross-correlation func- 
tions 

r,,(t) = {S^{to)S,{to + t)), (55) 

where i,j G {a.b} and Sa{t) is the plaquette operator creating a O"*" glueball from the 
vacuum, and Sb{t) is the operator creating a a-fg. If there is a non-zero mixing the 
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hermitean correlation matrix Fj^ would not be diagonal. More generally one defines 

where a; is a real valued parameter. Diagonalizing A{t) yields two eigenvalues, which are 
dominated by the lowest masses at large times t [38] : 

Ao(t) = /o(cj)e-™»*{l + 0(e-("i-™o)*)} (57) 
Ai(i) = /i(u;)e-"*i*{l + 0(e-^™i*)}, Ami = min(mi - mo, m2 - mi). (58) 

By tuning u the statistical errors can be minimized. The masses rriQ and mi belong to 
the two lightest physical states in this channel. The mixing angle 6(t) is defined to be the 
angle between the eigenvector Vo{t) corresponding to Aq and the vector (1,0). For large t 
one should observe a plateau where the mixing angle is constant and independent of cu. 
We have determined the mixing angle in the 0"*" channel from our Monte Carlo data. If 



uj takes its optimal value Uq = yTaa/^bb [36], the errors are smallest. Figure 13 shows the 
mixing ang le e{t) for this choice of u;. On the 8^ • 16 lattice for K = 0.19 and K = 0.1925 
as well as on the 12^ • 24 lattice for K — 0.1925 the result is consistent with zero within 
rather small errors. So there is no mixing between the glueball and the a-fo state. It 
might be possible that mixing only becomes visible in the close vicinity of the critical 
line corresponding to zero gluino mass, where supersymmetry is nearly restored. On the 
other hand, the effective action of [5] does not necessarily require a non-zero mixing to 
be present. 



5 Summary and outlook 

The numerical Monte Carlo simulations presented in [3] and this paper are the first 
calculations of this kind in a Yang-Mills theory with light gluinos. Therefore an essential 
part of our work had to be invested in algorithmic studies and parameter tuning. 

The two-step multi-bosonic algorithm, after appropriate tuning, turned out to be 
reliable and showed a satisfactory performance in the present case which is described 
by a fiavour number Nf — ^ oi fermions in the adjoint representation. We showed 
that the sign of the Pfaffian appearing in a path integral formulation of gluinos can be 
taken into account, but does not practically infiuence the results in the investigated range 
of parameters. Since the two-step multi-bosonic algorithm can also be applied for any 
number of fermion fiavours in the fundamental representation, an interesting physical 
application would be, for instance, QCD with three light fiavours of quarks. On the basis 
of our positive experience with the algorithm we expect that it would also work well in 
that case. 
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Concerning parameter tuning in the lattice action, the problem is to find a region of 
bare parameter space where the gluino is light and where the lattice spacing is appropriate 
for feasible numerical simulations. Our strategy was to start at the lower end of the 
approximate scaling region in pure SU(2) lattice gauge theory ai (3 — 2.3 and to increase 
the hopping parameter K as long as substantial effects of virtual dynamical gluino loops 
appear. It is expected that these effects decrease the lattice spacing due to the difference 
of the Callan-Symanzik /3-functions with and without light gluinos. The observed effect 
is mainly an overall renormalization of a. The change of dimensionless ratios of masses 
and string tension are only moderate up to < 0.1925, where most of our simulations 
were performed. 

Increasing K further is getting more difficult from the algorithmic point of view be- 
cause the smallest eigenvalues of the fermion matrix are becoming really small. In spite of 
this our algorithm still performed reasonably well. A search in the range up to X < 0.20 
revealed first evidence for a first order phase transition expected to occur at zero gluino 
mass [3]. Our present estimate for the location of this phase transition, at /3 = 2.3, is 
Kq = 0.1955(5). This gives for the bare gluino mass in lattice units 



1 



1 1 



(59) 



a value amo — 0.04 a,t K — 0.1925. With the value of the string tension in (46) we get 
fnol^fa ~ 0.2. Using QCD-units and neglecting the mass renormalization factor of 
the order of 1 this corresponds to a light gluino mass of about 100 MeT^. Of course, this 
can only serve as an order of magnitude guide because SYM and QCD are after all two 
different theories. In order to connect mo to, say, A^g^ one had to perform a calculation 
as in [39] with massless gluinos. 



Table 8: Masses of the light bound states at (3 — 2.3 in lattice units. 



lattice: L,T 


K 


0+ 


am% 


0+ 

am% 


am% 


am-yg 


6,12 


0.18 


0.95(10) 




1.49(13) 


1.155(11) 


1.93(5) 


6,12 


0.185 


0.85(6) 


1.5(3) 


1.11(17) 


0.941(8) 


1.39(8) 


8,16 


0.19 


0.75(6) 




1.20(22) 


0.725(20) 


0.87(13) 


8,16 


0.1925 


0.63(5) 


1.1(1) 


0.81(17) 


0.551(17) 


0.82(18) 


12,24 


0.1925 


0.53(10) 




1.00(13) 


0.48(5) 


0.93(8) 



Having an algorithm and knowing the interesting range of parameters in the lattice 
action one can start to perform numerical simulations for determining the spectrum of 
states and other physically interesting features. The properties of the lightest states are 
obviously quite interesting because the construction of low energy effective actions [4, 5] 
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is based on the assumptions about the relevant composite field variables. An important 
constraint on the spectrum is that in the limit of zero gluino mass, where supersymmetry 
is expected, the particle states should occur in supermultiplets with degenerate mass. A 
collection of our present results on the hghtest states is displayed in table 8 and figure 14. 

As one can see, for the lightest gluino masses (highest hopping parameters K) the 
bound state masses can be arranged into two groups. The lightest states are the 0~ 
gluinoball [a-r]') and the O"*" glucball. At K = 0.1925 these are in lattice units both near 
am ~ 0.5. The other group of states is at K = 0.1925 near am ~ 1.0 and consists of 
the 0"*" gluinoball, the 0~ glueball and the spin-| gluino-glueball. As shown by figure 13, 
there is practically no mixing between the 0^-states in the two groups. The disturbing 
fact concerning supersymmetry is that there is apparently no spin-| state in the lower 
mass group. We saw this problem already at early stages of our project and hence paid 
specific attention to a lighter spin-| gluino-glueball state, but we did not find it. 

There arc several possible explanations for this. Perhaps we are not yet close enough 
to the supersymmetric limit and therefore the spectrum does not yet look like a weakly 
broken supersymmetric spectrum. Another possibility is that we are missing the other 
spin- 1 state because our choice of interpolating fields is not appropriate. One can, for 
instance, think about spin-| gluinoballs made out of three gluinos which appear at strong 
coupling [40] and were not exploited in our simulations. Nevertheless, even if the spin-| 
state completing the lightest supermultiplet would be dominated by three gluinos, the 
emerging structure of the two light supermultiplets would be surprising. Finally one 
can think about possible finite volume effects and the effect of lattice artifacts breaking 
supersymmetry at finite lattice spacing. Without further numerical simulations we cannot 
exclude this last possibility but we believe that it is unlikely on basis of the experience in 
pure gauge theories. The product of the lattice spacing with the square-root of the string 
tension is at X = 0.1925 given by ay/a ~ 0.17. In pure SU(2) gauge theory [33] we have 
a similar value at /3 ~ 2.5 — 2.6 which is within the region of reasonably good scaling. 
As discussed in section 3, the spatial volume extension of our 12^ lattice at X = 0.1925 
is about 1 fm in QCD units. This is almost certainly not large enough and therefore 
there are important finite volume effects to be expected, but the qualitative features of 
the bound state spectrum should already be visible in such volumes. 

We leave this puzzle for further investigations. The most important outcome of this 
first Super- Yang-Mills simulation with light gluinos is that the numerical Monte Carlo cal- 
culations arc definitely possible with present-day techniques and can certainly contribute 
to the better understanding of the low energy non-perturbative dynamics of supersym- 
metric gauge theories. 
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Appendix 



A Least-squares optimized polynomials 

Least-squares optimization provides a general and flexible framework for obtaining the 
necessary optimized polynomials in multi-bosonic fermion algorithms. By exploiting dif- 
ferent weight functions this framework is well suited to fulfill rather different requirements. 

In the first part of this appendix the basic formulae from [17] are collected. In the 
second part a simple example is considered: in case of an appropriately chosen weight 
function the least-squares optimized polynomials for the approximation of the function 
x~°' are expressed in terms of Jacobi polynomials. 



A.l Definition and basic relations 

The general theory of least-squares optimized polynomial approximations can be inferred 
from the literature [41, 42]. Here we introduce the basic formulae in the way it has been 
done in [17] for the specific needs of multi-bosonic fermion algorithms. We shall keep the 
notations there, apart from a few changes which allow for more generality. 

We want to approximate the real function f{x) in the interval x e [e. A] by a polyno- 
mial Pn{x) of degree n. The aim is to minimize the deviation norm 

5„ = ^N-l j^xwixf [fix) - Pn{x)f^' . (60) 

Here w{x) is an arbitrary real weight function and the overall normalization factor N^^\ 
can be chosen by convenience, for instance, as 

N,^^ = j\xw{xff{xf . (61) 

A typical example of functions to be approximated is f\x) = x~"^ / P{x) with a > and 
some polynomial P{x). The interval is usually such that < e < A. For optimizing the 
relative deviation one takes a weight function w{x) = f{x)~^. 

It turns out useful to introduce orthogonal polynomials ^^{x) (// = 0, 1, 2, . . .) satis- 
fying 

J dxw{xf^^{x)^^{x) = S^^^^q^ . (62) 
and expand the polynomial P„ (x) in terms of them: 

Pr,ix)^J2dn.Mx) ■ (63) 

Besides the normalization factor let us also introduce, for later purposes, the integrals 
Pu and Su by 

r\ p\ r\ 

q^ = dx'w{xy^t,{xY , Pv = dx'w{x)^^t,{x)^x , Sj, = / dx'w{xyx'' . (64) 
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It can be easily shown that the expansion coefficients dnu minimizing 5„ are indepen- 
dent of n and are given by 

dn^ = d^^ — , (65) 

where 

K = J dxw{xff{x)^^{x) . (66) 

The minimal value of is 

n 

61^1- N-l Yl d.K . (67) 

The above orthogonal polynomials satisfy three-term recurrence relations which are 
very useful for numerical evaluation. The first two of them with /j, = 0,1 are given by 

^o{x) = 1 , $i(x) = X - — . (68) 

So 

The higher order polynomials $^(a;) for = 2, 3, . . . can be obtained from the recurrence 
relation 

$^+i(x) = (x + /5^)<l>^(x)+7^_i$^_i(a;) , (/x = 1, 2, . . .) , (69) 
where the recurrence coefficients are given by 

Using these relations on can set up a recursive scheme for the computation of the 
orthogonal polynomials in terms of the basic integrals defined in (64). Defining the 
polynomial coefficients /^j, {0 < u < /i) by 

M^) = E U^^'-" (71) 

the above recurrence relations imply the normalization convention 

/^o = l, (// = 0,1,2,...) , (72) 

and one can easily show that and satisfy 

9m = E fl^<yS2n-u , Pn = J2 ft^^ {S2f,+l-u + f^,lS2n-u) ■ (73) 

The coefficients themselves can be calculated from /n = — Si/sq and (69) which gives 

//U+1,2 — ftj.,2 + f^fiffi,! + 7m- 1 ; 
//i+1,3 — fiJ.,3 + /3/^//^,2 + Ifi-lffi-l,! , 

"I" 1 ~l~ T//— i>M— 2 ) 

= PiJ.fiJ.,iJ. + 1 ix-if ■ (74) 
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The polynomial and recurrence coefficients are recursively determined by (72)- (74). The 
expansion coefficients for the optimized polynomial Pn{x) can be obtained from (65) and 

b^. = E/m- ['dxw{xff{x)x^-^ . (75) 
A. 2 A simple example: Jacobi polynomials 

The approximation interval [e, A] can be transformed to some standard interval, say, [—1,1] 
by the linear mapping 

^- ^""x-e' ' ^=|(A-6) + i(A + 6). (76) 

A weight factor (1 + 4)^(1 — with p,a > —1 corresponds in the original interval to the 
weight factor 

Taking, for instance, p = 2a, a = this weight is similar to the one for relative deviation 
from the function f[x) = x~", which would be just x^". In fact, for e = these are exactly 
the same and for small e the difference is negligible. The advantage of considering the 
weight factor in (77) is that the corresponding orthogonal polynomials are simply related 
to the Jacobi polynomials [43, 44], namely 



Our normalization convention (72) implies that 



2x - A - e' 
A — e 



(78) 



aM -(X- eV^^^^'^^'v^ r(p + + l)r(a + + l)r(p + a + + 1) 



The coefficients of the orthogonal polynomials are now given by 

p-u; \ ( p \ r(p + // + i)r(p- 

v-uj ]\uj ] r(p + //-a; + l)r(p + (7 + 2// + l) 



w=0 

In particular, we have 



(p + (j + 2) ■ 



rx^ = 1 , Ar = -6 - (A - 6) . (81) 



The coefficients /3, 7 in the recurrence relation (69) can be derived from the known recur- 
rence relations of the Jacobi polynomials: 

2^ ^ ^^2(p + a + 2//)(p + (7 + 2// + 2) ' 

Ap,'^) ^ _(x _ ,^2 p{p + p){a + p){p + a + p) 

^'^"^ ^ ' {p + a + 2p-l){p + a + 2pf{p + a + 2p + l)- ^ > 
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In order to obtain the expansion coefficients of the least-squares optimized polynomials 
one has to perform the integrals in (75). As an example, let us consider the function 
f{x) — x""" when the necessary integrals can be expressed by hypergeometric functions: 



dx {x-eY{\-xYx^'-''- 



.p+,+i,^_,_„r(p + i)r(a + i) ^/^^ 1. 0.1 ^ 



^ (A - 6)^+-+^A^— " r(^^^^^^3) (^a - + i., a + 1; p + <7 + 2; 1 - - j . (83) 
Let us now consider, for simplicity, only the case e = 0, when we obtain 

= i D^A^+P+^+M-g ^l/^ + + + i)r(Q + /i)r(p - a + i)r(a + /i + 1) 

^ ^ ' r(p + a + 2/i+ l)r(a)r(p + (T-a + /i + 2) ' ^ ' 

Combined with (65) and (79) this leads to 

^ ( 1 x -M-g r(p + g + 2// + 2)V{a + p)r(p -a + l) 

^ ^ ' p!r(p + p + i)r(Q;)r(p + (7-a + /i + 2) ■ ^ > 

These formulae can be used, for instance, for fractional inversion. For the parameters 
p, (7 the natural choice in this case is p = 2a, a — Q which corresponds to the optimization 
of the relative deviation from the function f{x) — x~°'. As we have seen in section A.l, 
the optimized polynomials are the truncated expansions of x~°' in terms of the Jacobi 
polynomials p(^°='°) . The Gegenbauer polynomials proposed in [45] for fractional inversion 
correspond to a different choice, namely p — a — a — ^. This is because of the relation 

i (2Q;)i (n + a + 2 j 

Note that for the simple case a = 1 we have here the Chebyshev polynomials of second 
kind: Cl{x) = Unix). 

In our present application we have to consider a = |. For the first polynomial P^J^ we 
could, for instance, use the Gegenbauer polynomials corresponding to p(~4'~4). (For 
p(2,3,4) need, of course, the polynomials introduced in [15] which approximate more 
complicated functions.) A numerical comparison shows, however, that the least squares 
optimized polynomials minimizing the relative deviation in the interval [e. A] are better 
than the Gegenbauer polynomials (see fig. 15): both approximations are similar at the 
lower end of the interval but otherwise the deviations of the former are by a factor of five 
smaller. 

The special case a = | is interesting for the numerical evaluation of the zero mass 
lattice action proposed by Neuberger [46]. In this case, in order to obtain the least- 
squares optimized relative deviation with weight function w{x) — x, the function x'^ 
has to be expanded in the Jacobi polynomials P^^'^\ Note that this is different both 
from the Chebyshev and the Legendre expansions applied in [47]. The former would 
correspond to take P^~h~^\ the latter to P^°'^\ The corresponding weight functions 
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would be [x{X — x)] 2 and 1, respectively. As a consequence of the divergence of the 

weight factor at x — 0, the Chebyshev expansion is not appropriate for an approximation 

(-- --) 

in an interval with e = 0. This can be immediately seen from the divergence of ^' ^ 
at a = I in (85). 

The advantage of the Jacobi polynomials appearing in these examples is that they 
are analytically known. The more general least-squares optimized polynomials defined 
in the previous subsection can also be numerically expanded in terms of them. This is 
sometimes more comfortable than the entirely numerical approach. 



B Determining the smallest eigenvalues 

For finding the smallest eigenvalues and the corresponding eigenvectors of the squared 
hermitean fermion matrix = Q^Q we apply the algorithm of Kalkrcuter and Simma 
[48] . Some modifications and the optimization with detailed tests have been described in 
[13]. Here we give a short summary for the readers convenience. 

The smallest eigenvalue of a general hermitean matrix H can be found by minimizing 
the Ritz functional 

. ^ . (87) 

Here the notation defined in (31) is used, with z* denoting the complex conjugate vector 
of z and {xy) = {xly) = (yx), where / is the unit matrix. The gradient of the Ritz 
functional is obviously 

For the conjugate gradient procedure we can choose a starting search direction pi — —gH{z) 
and the iteration is defined by a new approximation to the eigenvector 

Zi+i = Zi + aiPi , i = 1, 2, . . . . (89) 

The factor ctj is chosen at the minimum of ij>h{z) in the search direction p^. One can show 
that 

2 {p:Hz,) 



0£i 



(90) 



{z*Hzi) - {ptHpi) - ^[{zlHzi) - {plHpi)f + 4 {plHzi){z*Hp,)] ' 

Let us note that taking the positive sign in front of the square root gives the maximum, 
instead of the minimum. The other sign can be used for finding the maximal eigenvalue 
instead of the minimal one. In the iteration relation (89) the conjugate search direction 
Pi+i can be chosen according to [48] 



Pi+i = gH{zi+i) + Pi 



Pi Zi+i-— r- 



(91) 
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For the factor /3j one can take, according to the Fletcher- Reeves prescription, with = 

A = %^ (92) 
or alternatively, according to the Polak-Ribiere prescription, 

g ^ {9i+l9i+l) - {9i9i + l) /gON 

It turned out that in case of our fermion matrices the Polak-Ribiere version is 25% to 40% 
more efficient than the Fletcher- Reeves version proposed in [48] . In naive implementations 
of this iterative procedure numerical problems may occur due to the increasing length of 
the vector Zj. Since the Ritz functional is scale invariant, this problem can be avoided by 
rescaling, typically every 25 steps, as 

Pi Pi \l {z*Zi) , 9i^ 9i \l {z*Zi) . (94) 



Several smallest eigenvalues might be determined by applying the above conjugate 
gradient iteration subsequently to the projection into the orthogonal subspaces defined 

by ^ ^ 

Hk = P^HP^ , P^v = v^Y: V, {vlv) , (A; = 2, 3, . . .) . (95) 

i=l 

Here Vi denote the previously found normalized eigenvectors. This naive procedure be- 
comes numerically instable after a few eigenvalues because of the numerical errors in the 
projectors P^. One can stabilize and speed up this sequential search if one embeds it in 
an iterative scheme [48]. If one is interested in the kmax smallest eigenvalues then, after 
finding some approximation to i'i,f2, . . . ,Vk^^^ in a sequential search, the k^ax ® k^ax 
matrix 

M,j = iv*Hvj) (96) 

is diagonalizcd. For reasonable values of k^ax this is a small problem and the resulting 
new eigenvalues and the corresponding eigenvectors 

v[=T. (97) 

are better than Vi. Here denotes the eigenvectors of the matrix M. After this interme- 
diate diagonalization the sequential search with conjugate gradient iterations is continued. 

After the restarting of the sequential search it takes some time until the search di- 
rections of the conjugate gradient iterations become again optimal. Therefore it is not 
good to insert an intermediate diagonalization too often, especially at later stages when 
the final precision is approached. In our project a good performance could be achieved 
if between the i-th and {i -\- l)-th intermediate diagonalization the sequential search was 
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performed with (5 + lOi) conjugate gradient iterations. The apphcation of the projectors 
Pj^ becomes, even for moderate values of k, quite expensive. Since Pj^ projects out ap- 
proximate eigenspaces of if, it is not necessary to apply it at every conjugate gradient 
iteration. Tests show that it suffices to perform the projection only in the intermediate 
diagonalizations and, say, after every 25-th conjugate gradient iterations. 

The optimization of the Kalkreuter-Simma algorithm pays off very well [13]. It turns 
out that the number of necessary conjugate gradient iterations per eigenvalue is getting 
smaller and smaller with increasing values of kmax- Another feature is that the last few of 
the kmax eigenvalues are slowly converging. As a consequence, for computing more than 
kmax = 16 eigenvalues, it is advantageous to run the algorithm with say, 5% larger 

than kmax and stop the iteration if the smallest kmax eigenvalues satisfy the stopping 
criterion. 

C High performance CH — h for LGT simulations 

When starting to develop the software for the DESY-Miinstcr Collaboration we decided 
to take care for the reusability and flexibility of the codc^. It is well known that object 
oriented design and programming (OOD, OOP) helps to fulfill these needs [49]. A widely 
spread prejudice against OOP is bad performance. But new techniques like expression 
templates [50, 51, 52] or temporary base class idiom [53] encouraged us to use an object 
oriented approach for the software development. There were serveral reasons, which led 
to the decision of using C-t--|- in our project. It is the only object oriented (00) language 
available for high performance computers and it is a high efficient 00 language. Message 
Passing (MPI) and multithreading libraries (POSIX threads) are also usable with C-|--|-. 
With the help of templates C-l-4- also supports generic programming [54]. This feature 
allows one, for instance, to write a code template for the whole lattice gauge theory (LGT) 
simulation without specifying the gauge group. By simply providing a gauge group class 
which describes the basic functionality of the desired gauge group (SU(2), SU(3) etc), 
the compiler is able to generate code. Generic programming made possible to port our 
Super- Yang-Mills simulation program from SU(2) to SU(3) in less than one week. The 
only thing we had to add was a high efficient SU(3) class which consisted of approximately 
800 lines (less than 2% of the project). 

By using these techniques the efficiency of the simulation code stands and falls with 
an efficient vector class. A problem that arises almost always while overloading numerical 
operators of a vector class is the generation of temporary objects. This problem is not 
^For more informations see http://pauli.uni-muenster.de/'' spander/susy/phd.c-|-+.ps 
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only limited to C++ but also to FortranQO. As a simple example let us consider 

a = f + y+i' . (98) 

ti 

^ V ' 

t2=i{+Z 

— * 

t\ is generated in a function, which means on the stack. If the function is left this tem- 
porary object has to be copied away from the stack using the so called copy constructor. 
That means the copy constructor is used to generate a temporary copy from a temporary 
object. Furthermore the compiler may generate hidden temporary vectors using the copy 
constructor. 

A popular method to avoid unnecessary copying is reference counting [55]. But due 
to the additional level of indirection reference counting is efficient only for large vectors 
and suited for typical vector length of dynamical fermion simulations. The basic ideas of 
two other solutions which are working fine for both, small and large vectors are 

• Temporary Base Class Idiom 

— introduces an own class Tmp Vector for temporary vectors, 

— Tmp Vector construct/destructed shallowly, 

— operator+ (Vector &) returns a Tmp Vector, 

— operator+ (Tmp Vector) is implemented as TmpVector+=Vector, 

— disadvantage: four times more operators have to be overloaded. 

• Expression Templates 

— avoids temporary objects in the first place by automatically transforming 
vector u,v,w; u=v+w; at compile time (more or less) into 

for (int i(0); i < u.lengthO; ++i) 
u[i]=v[i] + w[i] ; 

using template meta programming (or compile time programs). 

On one hand it is desirable to implement a class for handling gamma matrices. On 
the other hand it is obvious that the gamma matrix multiplication has to be done at 
compile time rather than at run time. Otherwise a fermion matrix multiplication would 
proceed at a snail's pace. Two techniques exist to achieve this. 

• Lazy Evaluation for 7^'0 

— delays computation until the result is needed. 

— processes expressions like x + IijM^ ^ single task. 
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• Expression 7^7,/ 

— forces the compiler to perform this multiphcation a compile time using template 
met a programming [56] 

To test the efficiency of different vector classes we used a Monte-Carlo simulation of 
the two dimensional cr-model. This is a worst case test for a vector class because the vector 
length is three. The administration overhead caused by the introduction of a class can 
be huge compared to the performed operation. Generally the difference between the class 
libraries are small for larger vector sizes. As one can see in tabular (9) Blitz-|— |- (Expression 
templates) and NumArray.h (temporary base class) reach comparable speed to a hand- 
optimized Fortran?? implementation on a T3E-512/600. MV-|— I- uses reference counting 
which is not suitable for small vector sizes. The only commercial library math.h-l— |- 
surprisingly is the slowest. 



Fortran?? 


Blitz-F-F 


NumArray.h 


MV++ 


math.h-|— 1- 


4.81s 


4.93s 


5.21s 


40.1s 


69.1s 



Table 9: Runtimes of various vector classes for the simulation of the 2d 0(3) sym- 
metric non-linear (7 -model on the T3E-512/600 with the Cray C-h-h compiler. 

Usually larger object oriented programs break up into packages which are only loosely 
connected. Unfortunately C++ does not support the decomposition in modules as for 
example JAVA does. The package structure of the simulation is shown in diagram 16. 
The main ingredients are algorithms which act on fields. They make up 90% of the code. 
With the iterator pattern [5?] and an abstract I/O concept the corresponding code is 
hardware independent. It is suited for massive parallel, symmetric multiprocessor and for 
single CPU architectures. The hardware dependent objects, like iterators or I/O streams, 
should not be created by objects of these packages, but by the central object factory. 

Depending on the hardware on which the program is running, the object factory 
generates the suitable objects. The only hardware dependent components are the iterators 
and the I/O system. The question might arise why a dedicated I/O system for SMP is 
missing. The answer is that it is not needed. Only the algorithms really use more than 
one thread. When the algorithm is completed all threads are joined to a single one and 
this one uses native I/O routines. 
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Figures 
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Figure 1: The distribution of the correction factors in nine independent (parallel) 
sequences of configurations on 12^-24 lattice at (3 = 2.3, K = 0.1925. The considered 
configurations are separated by 50 updating cycles. The upper part shows the distri- 
bution and a Gaussian fit. In the lower part the independent lattices are separated 
by vertical lines. 
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Figure 2: The measurement correction for the a-pion propagator at zero distance. 
The exceptional configurations with small eigenvalues contribute strongly to the raw 
data. After correction these contributions are still important but of normal size. 
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Figure 3: The spectral flow of the hermitean fermion matrix Q for some specific 
configurations on 6^ • 12 lattice at {3 = 2.3. The value of K in the simulation is 
displayed by a vertical line. 
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Figure 4: The spectral flow of the hermitean fermion matrix Q for some configura- 
tions separated by 50 updating cycles on 6^ • 12 lattice at (5 — 2.3; K — 0.2. 
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Figure 5: The distribution of the smallest eigenvalues of the squared preconditioned 
fermion matrix Q'^ versus the non-preconditioned one on a 6^ ■ 12 lattice at 
(3 = 2.3, K = 0.196. 
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Figure 6: Relative deviation of the successive polynomial approximations of x^^^^ in 
the range of eigenvalues corresponding to our simulations in the 12^ • 24 lattice at 
K = 0.1925. (For parameters see table 2.) 
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Figure 7: Autocorrelation function and exponential fit for the gluino-glue propagator 
on one of the 12^ ■ 24 lattices run in parallel at [3 = 2.3, K = 0.1925. 
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Figure 8: Typical autocorrelation of the plaquette, with the exponential fit. The right 
graph shows the same data in a smaller interval. 
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Figure 9: Potential V{R, T) for 2, as a function of T on a 12^ • 24 lattice. The 
line is an exponential fit to the large T behaviour, fitted over the range 1 < T < 6. 



Lattice: S-'xie, K = 0.19 




Figure 10: The static quark potential V{R) on a 8'^ ■ 16 lattice at K = 0.19. The 
line is a fit with a Coulomb plus a linear term, fitted over the range 1 < R < 4. 
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Lattice: 12^x24, K = 0.1925 




Figure 12: Closed loop, which has been used to build the pseudoscalar glueball oper- 
ator. 
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Lattice: 12^x24, K^O.1925, (3 = 2.5 
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Figure 13: The mixing angle 6{t) in the 0+ channel on a 12^-24 lattice at K = 0.1925. 
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Figure 14: The lightest bound state masses in lattice units as function of the bare 
gluino mass parameter 1/K. The shaded area at K = 0.1955(5) is where zero gluino 
mass and supersymmetry are expected. 
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Figure 15: Comparing the polynomial approximations of x in the interval 
[e,A] = [0.003,3.7] at the order Ui = 32. P^^^{x) x^^l^ is shown for the least squares 
optimized polynomial minimizing the relative deviation (full line) and for the frac- 
tional inversion defined by the Gegenhauer polynomials with index a = ^ (dashed 
line). The interval is shown in three parts in order to display better the details. 
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Figure 16: UML packet structure diagram, showing the hardware dependent and 
independent parts of the project. 
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